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Final  Progress  Report 

Nonlinear  Solvers  for  Subsurface  Flow  Problems 

ARO  Grant  #  DAAD 19-99- 1-0 186 
May  1,  1999 -April  30,  2002 

1.  Problem  Studied.  This  project  was  for  basic  research  in  numerical  methods  for  solving 
linear  and  nonlinear  equations  and  related  problems  in  time-dependent  simulations.  The  motivating 
application  was  simulation  of  flow  through  porous  media,  in  particular  the  unsaturated  zone  in  the 
subsurface.  The  nature  of  the  nonlinearities  and  the  physical  properties  of  the  equations  make 
standard  methods  for  solution  of  equations  perform  in  ways  not  predicted  by  current  theory  [17]. 
Our  objective  was  to  understand  these  effects  and  design  and  analyze  new  solvers.  By  doing  this 
the  new  solvers  can  be  incorporated  in  production  codes  and  adapt  to  changes  in  the  solution  as  the 
simulation  progresses.  This  ability  to  adapt  is  crucial  to  an  efficient  simulation. 

The  basic  research  issues  in  nonlinear  equation  solvers  arise  from  the  nature  of  the  nonlineari¬ 
ties.  These  do  not  have  the  smoothness  properties  that  a  traditional  Newton’s  method  code  requires 
to  converge  in  the  usual  way.  Our  group  performed  research  directed  toward  understanding  of  how 
the  nonlinearity  can  be  approximated  in  a  way  that  improves  the  performance  of  the  solver  and, 
at  the  same  time,  does  not  affect  the  accuracy  of  the  simulation.  The  nonlinearities  also  affect 
the  way  in  which  the  nonlinear  solver  and  the  time -dependent  part  of  the  simulator  communicate. 
This  communication  is  important  for  temporal  adaption.  This  part  of  the  project  is  the  central 
component  in  Kathleen  Kavanagh’s  Ph.  D.  work. 

The  research  questions  for  linear  solvers  arise  from  the  size  of  the  problem.  In  general  terms, 
the  linear  solver  is  used  to  compute  an  approximation  to  the  Newton  step  for  the  nonlinear  solver. 
This  part  of  the  overall  simulation  takes  most  of  the  computational  time.  Only  a  large  distributed- 
memory  parallel  computer  can  store  the  data  for  large-scale  ground  and  surface-  water  simulations 
in  three  space  dimensions.  Linear  solvers  that  are  efficient  in  this  environment  and  can  effectively 
use  the  complex  data  structures  needed  for  an  adaptive  spatial  discretization  must  be  designed. 
In  collaboration  with  the  group  at  WES  we  have  designed  and  implemented  in  ADH  an  efficient 
iterative  method  for  large  linear  systems  [4,6-8].  This  method  has  given  good  parallel  scalability 
and  is  now  [8]  supported  by  theoretical  analysis. 

This  research  was  done  in  collaboration  with  a  group  (Charlie  Berger,  Stacy  Howington, 
and  Jackie  Hallberg)  at  the  US  Army  Engineer  Research  and  Development  Center  (ERDC).  The 
Adaptive  Hydrology  model  (ADH),  a  production  groundwater  modeling  code,  was  used  as  both  a 
testbed  for  the  algorithmic  work  in  the  project  and  as  a  source  for  new  research  topics.  The  PI  and 
his  students  have  an  very  productive  collaboration  A  DH  team. 

In  addition  to  the  group  at  ERDC,  the  PI  collaborates  a  group  led  by  Professor  C.  T.  Miller  at 
the  University  of  North  Carolina  on  work  related  to  this  project. 

2.  Most  Important  Results.  The  most  significant  results  from  the  project  are  connected  with 
the  development,  implementation,  analysis,  and  tech  transfer  of  a  two-level  Schwarz  precondi¬ 
tioner  for  Richards’  equation  (RE)  and  related  problems.  Richards’  equation  is  a  simple  model  of 
flow  in  the  unsaturated  zone.  In  three  space  dimensions,  letting  z  be  the  vertical  direction  and  V 
the  spatial  gradient  operator,  the  pressure  head  form  of  RE  is 

(2.1)  [c(V0  +  SA(v)]  ^  =  V  •  (/iW’)Ve)  + 
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In  (2.1),  xp  is  pressure  head,  r(t;)  =  dO /dp'  is  the  specific  moisture  capacity,  9(p>)  is  the  volumetric 
fraction  of  water,  Ss  is  the  specific  storage,  Sa(p>)  =  9(p>)/n  is  the  aqueous-phase  saturation,  n  is 
the  porosity,  and  Kpw)  is  the  hydraulic  conductivity.  We  will  assume  that  appropriate  initial  and 
boundary  conditions  have  been  imposed. 

he  equation  must  be  closed  with  constitutive  equations  for  9pw)  and  K(p>).  One  common  way 
to  do  this  is  to  define  the  effective  saturation  Se  with  the  van  Geneuchten  formula  [19], 


(2.2) 


Se(p’) 


6-6r 
9S  —  9r 


(1  I  n, m"  ,  p’  <  0 
1,  >  0  ' 


In  (2.2),  9r  is  the  residual  volumetric  water  content,  9S  is  the  saturated  volumetric  water  content, 
nu  is  an  experimentally  determined  measure  of  pore  size  uniformity,  mv  =  1  —  1  jnv,  and  av  is 
an  experimentally-determined  coefficient  that  is  related  to  the  mean  pore  size.  Note  that  (2.2)  also 
defines  9  as  a  function  of  p>.  If  the  van  Geneuchten  formula  is  used  for  the  saturation,  it  is  standard 
to  define  the  conductivity  with  the  Mualem  model,  [15], 

(2.3)  K(4>)  =  KsSlJ 2  [l  -  (l  -  S1Jm")mv]2 

where  Ks  is  the  water-saturated  hydraulic  conductivity.  One  can  see  from  (2.2)  and  (2.3)  that  the 
nonlinearity  is  not  differentiable  if  1  <  n„  <  2.  Values  of  nv  in  this  range  are  physically  realistic 
and  solvers  must  be  prepared  to  deal  with  nonsmooth  nonlinearities. 

Discretization  in  time  and  space  of  RE  leads  to  a  sequence  of  nonlinear  equations  that  must 
be  solved  at  each  time  step.  To  see  this  in  a  simple  way  we  define  the  accumulation  term 


A(ip)  =  c(p’)  +  S8Sa(p>), 


and  suppress  the  spatial  variables  to  obtain 

A{P’)^r  =  AT  {P’) 

where  the  nonlinear  term  N  contains  all  spatial  derivatives.  The  semi-discretized  problem  (i.  e. 
discrete  in  space  only)  has  the  same  form, 

(2.4)  A(u)ut  =  J\f{u), 

where  N  is  the  discrete  form  of  J\f  and  u  the  approximation  to  pj.  Our  work  on  this  project  seeks 
efficient  solvers  for  (2.4). 

One  must  integrate  (2.4)  implicitly  in  time.  If,  for  example  one  uses  the  backward  Euler 
method  one  obtains 

(2.5)  A{un+1){un+1  -  un )  =  hnN(un+1) 

which  leads  to  a  nonlinear  equation  that  must  be  solved  to  advance  in  time 

(2.6)  F(u)  =  A(u)(u  -  un)  -  hnN(u)  =  0. 


We  solved  (2.6)  with  a  Newton-iterative  method.  This  means  that  we  approximate  the  Newton 

step 

S  =  -Fvr^K) 
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F\uc)s  =  - F(uc ), 


by  solving  the  linear  equation, 

(2.7) 

with  an  iterative  method. 

The  problems  of  interest  here  are  too  large  for  a  Jacobian  factorization  to  be  computed  and 
stored.  Hence  (2.7)  must  be  solved  by  an  linear  (inner)  iterative  method.  The  termination  criteria 
for  the  iterative  method  is  typically  coupled  to  the  state  of  the  nonlinear  (outer)  iteration  by  the 
inexact  Newton  condition  [1, 13], 

(2.8)  ||771,(ttc)s  +  F(uc)\\  <  ryc||F(ttc)||. 

This  is  simply  the  standard  relative  linear  residual  termination  criterion  for  iterative  methods. 

In  our  computational  work  on  RE  [4, 6,7]  and  two-phase  flow  [12]  we  use  a  preconditioned 
Krylov  method  for  the  linear  equation  and  terminate  that  inner  iteration  when  (2.8)  holds.  The 
problems  are  large  enough  that  distributed  memory  computers  are  needed  both  to  store  the  data  and 
to  do  the  computations  rapidly.  The  Jacobian  matrix  is  nonsymmetric  for  most  of  the  problems  we 
consider  and  we  solve  the  nonsymmetric  linear  equation  for  the  Newton  step  with  a  preconditioned 
Bi-CGSTAB  [18]  linear  iteration. 

Let  Ax  =  b  be  the  linear  system  to  be  solved.  We  seek  a  (left)  preconditioner  M  such  that  the 
condition  number 

k{MA)  =  ||MA||||(MA)-1|| 

is  significantly  smaller  than  that  of  A.  Krylov  methods  usually  perform  better  if  the  condition 
of  the  linear  system  can  be  improved.  The  preconditioner  is  used  by  applying  the  method  to  the 
preconditioned  system  M  Ax  =  Mb.  Krylov  methods  take  one  or  two  matrix- vector  products  per 
iteration.  The  matrix-vector  product  is  usually  the  most  significant  cost  in  cpu  time.  Hence  a 
good  preconditioner,  by  reducing  the  number  of  linear  iterations,  will  reduce  the  cost  of  the  solve. 
However,  one  must  keep  in  mind  that  the  M A- vector  product  could  be  much  more  expensive 
than  the  1- vector  product.  Hence,  the  preconditioner  is  useful  if  the  reduction  in  the  number  of 
iterations  offsets  the  increased  cost  from  the  application  of  the  preconditioner  with  each  matrix- 
vector  product. 

The  combination  of  a  Newton-Krylov  method  with  a  Schwarz  domain  decomposition  pre¬ 
conditioner  is  called  a  Newton-Krylov-Schwarz  (NKS)  method  [14].  We  have  implemented  both 
one  and  two  level  Schwarz  preconditioners  in  ADH.  Both  of  these  preconditioners  are  domain  de¬ 
composition  preconditioners,  which  means  that  the  original  physical  domain  is  split  into  several 
subdomains,  and  the  solutions  of  the  original  problem  restricted  to  the  subdomains  are  combined 
to  form  the  preconditioner  for  the  original  system. 

We  now  describe  additive  the  kind  of  Schwarz  preconditioners  we  use  in  ADH.  If  we  define 
a  matrix  R,  to  be  the  restriction  matrix  for  subdomain  i  so  that  II,  =  [0  I  0],  where  /  is  an 
n,;  x  n,  identity  matrix  and  n,  is  the  size  of  subdomain  i,  then  the  one-level  additive  Schwarz 
preconditioner  can  be  written  as 


M  =  t  ( RiARj)~ 'Ri 

2  =  1 


where  p  is  the  number  of  subdomains. 
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The  coarse  mesh  component  of  the  preconditioner  is  formed  by  defining  one  aggregate  ele¬ 
ment  per  subdomain.  The  resulting  coarse  mesh  basis  function  is  constant  except  in  the  elements 
shared  between  subdomains.  The  contribution  of  the  subdomain  to  the  coarse  matrix  is  computed 
locally  and  then  communicated  to  all  of  the  processors.  Thus  every  processor  solves  the  coarse 
mesh  problem.  The  subdomain  solves  are  performed  using  a  profile  solver  [2]  and  the  coarse  grid 
problem  is  solved  using  a  dense  LU  factorization.  The  two-level  additive  Schwarz  preconditioner 
is  formed  by  adding  the  coarse  mesh  problem  to  the  one-level  preconditioner,  so  that 

M  =  Rl  (R0ARl)~L  i?0  +  E  RJ  (RiARjy1  R, 

i 

where  R0  and  Rq  are  the  restriction  and  interpolation  operators  from  the  fine  to  coarse  meshes, 
respectively. 

If  h  is  the  scale  of  the  fine  mesh  and  H  the  scale  of  the  subdomains,  we  have  shown,  at  least 
for  model  elliptic  problems  [5,8],  that  the  two-level  preconditioner  satisfies, 

(2.9)  k(MA)  <Cuj{  1  +  ( H/hf ). 

In  practical  terms,  (2.9)  means  that  if  the  number  of  subdomains  grows  as  the  mesh  is  refined, 
then  the  condition  number  of  M  A  will  be  independent  of  h,  hence  the  number  of  iterations  needed 
to  satisfy  (2.8)  should  also  be  independent  of  h  and  the  scalability  of  the  solver  should  be  good. 
While  the  link  between  condition  number  and  iteration  count  is  only  a  heuristic  for  nonsymmetric 
systems  [3, 13, 16],  the  computational  results  in  [4, 6-8, 12]  show  good-to-excellent  scalability  in 
terms  of  iteration  counts  and  nonlinear  function  evaluations. 

To  illustrate  this  scalability,  we  present  a  table  of  iteration  statistics  from  [8].  The  findings 
in  [6,  12]  are  similar,  with  the  one-level  method  performing  less  well  in  the  three  dimensional 
simulations  from  [6].  The  data  are  from  a  temporal  integration  of  RE  in  two  space  dimensions, 
using  a  regular  spatial  mesh.  In  the  table  we  report  on  the  average  number  of  linear  iterations  per 
Newton  step  over  a  complete  simulation.  The  linear  solver  was  BiCGSTAB  [18]  and  the  temporal 
integration  was  done  with  the  BDF  code  from  [10].  Table  2.1  shows  how  the  performance  of 
the  preconditioned  linear  solver  depends  on  H  and  h.  The  constant  numbers  along  the  diagonals 
indicate  almost  perfect  scalability. 


Table  2.1 

RE  Iteration  Statistics,  2-level  Schwarz  Preconditioner 


In  [12]  we  apply  the  solver  technology  developed  for  RE  in  [5,  6,  8]  to  a  model  [9, 11]  for 
two  phase  flow.  The  model  is  a  coupled  system  of  equations  of  the  same  form  as  (2.1).  The 
preconditioning  and  temporal  integration  methods  that  have  been  designed  for  Richards’  equation 
also  performed  will  in  this  context. 
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In  [4,7]  the  group  at  ERDC  and  the  Pi’s  group  applied  the  multilevel  preconditioner  to  prob¬ 
lems  in  surface  water  and  in  ground-surface  water  interaction.  The  long-term  goal  is  that  the 
preconditioner  be  applicable  to  the  entire  range  of  problems  that  ADH  is  intended  to  solve. 
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W.  G.  Gray,  and  G.  F.  Pinder,  eds.,  Amsterdam,  2002,  Elsevier,  pp.  249-256. 
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Krylov-Schwarz  method  for  hydrology,  in  Parallel  Computational  Fluid  Dynamics  1999, 
D.  E.  Keyes,  A.  Ecer,  J.  Periaux,  and  N.  Satofuka,  eds.,  North  Holland,  2000,  pp.  257-264. 
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Carolina,  2000. 

4.  Scientific  Personnel. 

1.  C.  T.  Kelley:  Principal  Investigator 

2.  E.  W.  Jenkins:  Graduate  Student 

Ph.  D.  in  Mathematics,  August,  2000 

3.  K.  R.  Kavanagh:  Graduate  Student 

MS  in  Mathematics,  December,  2000;  Ph.  D.  expected,  2003 
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